clc;
clear all;
close all;

filename = 'Scan_3.jpg';

im = imread(filename);
im = im2double(im);
bw = rgb2gray(im);

%% resize
scale = 0.10;
bw = imresize(bw, scale);
bw = imadjust(bw);

[imH, imW] = size(bw);
imSz = imH*imW;

%% histogram
h = hist(bw(:), 256);
% bar(h);

barrier = 200;
sigma = 10/256;
[~, hist_max] = max(h(barrier:end));
hist_max = (barrier+hist_max-1)/256;

%% threshold
bw(bw>(hist_max-sigma)) = 0;
bw(bw~=0)=1;

%% erode
se = strel('disk',round(scale*15));
bw = imerode(bw,se);

%% filter out weak connection
conn_max_size = 30*scale;  % 20 pixels
h = fspecial('disk', round(conn_max_size*3));
blurred = imfilter(bw,h,0);
bw(blurred<0.4) = 0;

%% segmentation
seg = bwconncomp(bw, 4);
labeled = labelmatrix(seg);
RGB_label = label2rgb(labeled, @spring, 'c', 'shuffle');

%% find big components
lens = zeros(1, length(seg.PixelIdxList));
for ii=1:length(seg.PixelIdxList)
    lens(ii) = length(seg.PixelIdxList{ii});
end;

med_size = median(lens)*imSz*0.01;
pixIndList = {};
for ii=1:length(seg.PixelIdxList)
    len = length(seg.PixelIdxList{ii});
    if (len > med_size)
        pixIndList = [pixIndList, seg.PixelIdxList{ii}];
    end;
end;

imshow(im); hold on;
%% separate and send for inspection
for ii=1:length(pixIndList)
    inds = pixIndList{ii};
    im2 = bw*0;
    im2(inds) = 1;
    im2 = imdilate(im2,se);
    P = FindQuadrature(im2);
    plot(1/scale*P(1,:), 1/scale*P(2,:), 'r');
    cropped_im = CropImage (im, 1/scale*P);
end;
hold off;


